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Abstract 

The  most  common  theoretical  approach  to  model 
the  interactions  in  a  biochemical  process  is  through 
chemical  reactions.  Often  for  these  reactions,  the 
dynamics  of  the  first  M-order  statistical  moments 
of  the  species  populations  do  not  form  a  closed 
system  of  differential  equations,  in  the  sense  that  the 
time-derivatives  of  first  M-order  moments  generally 
depend  on  moments  of  order  higher  than  M.  However, 
for  analysis  purposes,  these  dynamics  are  often  made 
to  be  closed  by  approximating  the  needed  derivatives 
of  the  first  M-order  moments  by  nonlinear  functions 
of  the  same  moments.  These  functions  are  called  the 
moment  closure  functions. 

This  paper  presents  a  systematic  procedure  to 
construct  these  moment  closure  functions.  This  is  done 
by  first  assuming  that  they  exhibit  a  certain  separable 
form,  and  then  matching  time  derivatives  of  the  exact 
(not  closed)  moment  equations  with  that  of  the  approx¬ 
imate  (closed)  equations  for  some  initial  time  and  set 
of  initial  conditions.  Using  these  results  a  stochastic 
model  for  gene  expression  is  investigated.  We  show 
that  in  gene  expression  mechanisms,  in  which  a  protein 
inhibits  its  own  transcription,  the  resulting  negative 
feedback  reduces  stochastic  variations  in  the  protein 
populations. 

1  Introduction 

This  paper  presents  a  systematic  procedure  for  con¬ 
structing  approximate  stochastic  models  for  chemical 
reactions  used  for  modeling  biochemical  processes 
such  as  gene  regulatory  networks.  Such  models  are 
motivated  by  the  recent  work  [1,  2]  which  provide 
considerable  experimental  evidence  for  stochastic 
fluctuations  in  gene  expression  and  regulation  and  may 


account  for  the  large  amounts  of  cell  to  cell  variation 
observed  in  genetically  identical  cells  exposed  to  the 
same  environment  conditions  [3,  4].  Furthermore, 
studies  of  engineered  genetic  circuits  designed  to 
act  as  toggle  switches  or  oscillators  have  revealed 
large  stochastic  effects.  Stochastically  is  therefore  an 
inherent  feature  of  biological  dynamics  and  developing 
stochastic  models  which  capture  this  stochastically 
have  become  increasingly  important.  Analysis  of  such 
model  not  only  helps  us  to  discover  the  benefits  that 
biology  draws  from  stochasticity  but  also  helps  to 
explore  its  use  for  other  applications,  for  example, 
design  of  stochastic  bio-inspired  decision  algorithms 
to  control  the  motion  of  networks  of  artificial  mo¬ 
bile  agents  such  as  small  autonomous  micro-UAVs 
or  UGVs  involved  in  surveillance  and/or  tracking 
missions.  Although  our  work  focuses  on  biochemical 
reactions,  the  modeling  tools  developed  in  this  paper 
can  be  applied  to  a  very  general  class  of  stochastic 
systems,  in  particular  Markovian  systems  whose 
state  evolves  due  to  events  triggered  by  stochastic 
processes.  Other  examples  of  such  systems  include 
ecological  systems  where  these  events  correspond  to 
births,  recruitment  or  deaths  and  in  networked  control 
systems  where  the  events  correspond  to  exchanges  of 
messages,  communication  faults,  etc. 

The  most  common  theoretical  approach  is  to  model 
the  interactions  in  a  biochemical  process  as  chemical 
reactions.  The  stochasticity  in  the  biochemical  process 
is  then  captured  using  the  stochastic  formulation  of 
chemical  kinetics  which  treats  the  reactions  as  proba¬ 
bilistic  events.  The  time  evolution  of  the  system  is  then 
described  by  a  single  equation  for  a  probability  density 
function,  where  time  and  species  populations  appear 
as  independent  variables,  called  the  Chemical  Master 
Equation  (CME)  [5,  6].  This  equation  can  only  be 
solved  analytically  for  relatively  few,  highly  idealized 
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cases  and  generally  Monte  Carlo  simulation  tech¬ 
niques  are  used  which  are  also  known  as  the  Stochastic 
Simulation  Algorithm  (SSA)  to  study  stochasticity 
in  bio-chemical  reactions  [7,  8,  4,  9,  10].  Since  one 
is  often  interested  in  only  the  first  and  second  order 
statistical  moments  for  the  number  of  molecules  of 
the  different  species  involved,  much  effort  can  be 
saved  by  applying  approximate  methods  to  produce 
these  low-order  moments,  without  actually  having  to 
solve  for  the  probability  density  function.  Various 
such  approximate  methods  have  been  developed,  for 
example,  using  the  Fokker-Plank  approximation,  ex¬ 
panding  the  Master  equation,  etc  [5,  11].  In  this  paper, 
an  alternative  approximate  method  for  estimating 
lower-order  moments  is  introduced  using  moment 
closure  techniques. 

We  show  that  the  biochemical  reactions  can  be 
conveniently  modeled  as  a  Stochastic  Hybrid  System 
(SHS),  the  state  of  which  are  the  populations  of 
different  species  involved  in  the  reactions  [12],  Then, 
the  time  evolution  of  the  moments  of  the  population 
is  obtained  using  results  from  the  SHS  literature.  It 
has  been  shown  that  for  reactions  with  more  than  one 
reactants,  time  evolution  of  the  first  M  order  moments 
of  the  population  is  not  closed,  in  the  sense  that  it 
depends  on  moments  of  order  higher  than  M.  For 
analysis  purpose,  the  time  evolution  of  the  first  M 
order  moments  is  made  to  be  closed  by  approximating 
these  higher  order  moments  as  a  nonlinear  function 
of  moments  up  to  order  M,  which  we  refer  to  as 
the  moment  closure  function.  This  paper  presents 
explicit  formulas  to  construct  these  moment  closure 
functions  using  the  recently  introduced  techniques 
[13,  14]  of  matching  time  derivatives  between  the  the 
exact  (not  closed)  moment  equations  with  that  of  the 
approximate  (closed)  equations  for  some  initial  time 
and  set  of  initial  conditions.  The  striking  feature  of 
these  formulas  are  that  they  are  independent  of  the 
reaction  parameters  (reaction  rates  and  stoichiometry), 
moreover,  the  accuracy  of  the  approximation  can  be 
improved  by  increasing  M. 

These  closed  moment  equations  provide  time  evolution 
of  lower  order  moments  for  populations  of  species 
involved  in  the  biochemical  reactions.  Apart  from 
providing  fast  simulation  times  and  lesser  computation 
burden  compared  to  Monte  Carlo  simulations  these 
approximate  models  also  open  the  doors  to  other  types 
of  analysis  tools,  for  example,  sensitivity  analysis 
of  chemical  master  equation  with  respect  to  reaction 


parameters.  However,  they  provide  lesser  information 
about  the  probability  distribution  as  compared  to 
Monte  Carlo  simulations.  To  illustrate  the  applicability 
of  our  results,  we  investigate  a  stochastic  model  for 
gene  expression.  We  show  that  a  negative  feedback 
mechanism  caused  by  the  expressed  protein  inhibiting 
its  own  transcription,  reduces  stochastic  variations  in 
the  protein  populations. 

2  Stochastic  Modeling  of  Chemi¬ 
cally  Reacting  Systems 

Consider  a  system  of  n  species  Xj,  j  £  { 1 , . . . ,  n}  inside 
a  fixed  volume  V  involved  in  K  reactions  of  the  form 

Rj  :  M;  i  X |  +  . . .  +  Ui„Xn  — >  vnX\  +  . . .  +  vmXn  +  * 

(1) 

for  all  i  £  {1, . . . ,  AT},  where  u/j  £  N>o  is  the  stoichiom¬ 
etry  associated  with  the  jth  reactant  in  the  ith  reaction 
and  Vfj  £  N>o  is  the  stoichiometry  associated  with  the 
jth  product  in  the  ith  reaction,  and  *  represents  products 
other  than  the  species  Xj.  As  all  chemical  reactions  oc¬ 
cur  in  a  series  of  elementary  reactions  [15],  which  are 
generally  uni-  or  bi-molecular,  we  assume 

Uji  +  . . .  +  Ui„  <  2,  Vi  £  {1, . . . , K},  (2) 

and  hence,  we  only  allow  reactions  which  have  the 
form  given  in  the  first  column  of  Table  1 .  The  reac¬ 
tion  parameter  c;  characterizes  the  reaction  R,  and,  to¬ 
gether  with  the  stoichiometry,  defines  the  probability 
that  a  particular  reaction  takes  place  in  an  “infinitesi¬ 
mal”  time  interval  (t,t  +  dt\.  This  probability  is  given 
by  the  product  Cjhjdt  where  hj  is  the  number  of  distinct 
molecular  reactant  combinations  present  in  V  at  time  t 
for  the  reaction  R,  and  c,dt  is  the  probability  that  a  par¬ 
ticular  reactant  combination  of  R\  will  actually  react  on 
(t,t+dt\.  The  number,  h,  depends  both  on  the  reac¬ 
tants  stoichiometry  Ujj  in  R,  and  on  the  number  of  reac¬ 
tant  molecules  in  V.  Table  1  shows  the  value  of  h,  for 
different  reaction  types  [7],  In  this  table  and  in  the  se¬ 
quel,  we  denote  by  x;,  the  number  of  molecules  of  the 
species  Xj  in  the  volume  V  and  x  :=  [xi , . . .  ,x„]r  £  R". 

To  model  the  time  evolution  of  the  number  of 
molecules  xi,  X2,  . . . ,  x„,  a  special  class  of  Stochastic 
Hybrid  Systems  (SHS)  were  introduced  in  [12]  .  More 
specifically,  to  fit  the  framework  of  our  problem,  these 
system  are  characterized  by  trivial  dynamics 

x  =  0,  x=  [xi,...,x„]r,  (3) 
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Table  1:  h,(x)  for  different  reaction  types. 


Reaction  R; 

hi(x) 

Xj  >  * 

xj 

Xj+X i— ►*,  (tfj) 

XjX, 

2 Xj  —  * 

3X/(X/  —  !) 

a  family  of  K  reset  maps 

x  =  fc(x"),  (4) 

and  a  corresponding  family  of  K  transition  intensities 

A ,(x),  A,- :  R”— » [0,<=o)  (5) 

for  all  /  G  { 1 , . . . , /T}.  Each  of  the  reset  maps  0;(x),  and 
corresponding  transition  intensities  A,(x)  are  uniquely 
defined  by  the  ith  reaction  and  given  by 


If  any  of  the  reactions  have  more  than  one  reactants,  i.e, 
has  the  form  Xj  +  Xt  *  for  which  A,  (x)  is  quadratic 

in  x,  then,  from  (9)  it  can  been  shown  that  the  time 
derivative  of  a  moment  of  order  m* ,  is  given  by  a  linear 
combination  of  moments  of  orders  upto  m*  +  1  [17], 
Hence,  if  one  stacks  all  moments  in  an  infinite  vector 

M~  =  [M(ml),M(m2),---]T,  m PGN">0l  Vp G  {1,2,...} 

(10) 

its  dynamics  can  be  written  as 

[loo  =  Aoofloo  ,  (11) 

for  some  infinite  matrix  A„.  As  the  above  infinite  di¬ 
mensional  system  cannot  be  solved  analytically,  we 
truncate  (11)  by  creating  a  vector 

p  =  [ju(ml),fi(m2),-  ■  •  ,p(mk)]T  G  Rk  (12) 


X  0,-(x) 


xt  —  M/1  +  V/l 
x2  —  M/2  +  V/2 


A/(x)  =  Cihj(x)  (6) 


Min  H-  Vin 


for  all  i  G  {1, ...  ,K}.  In  essence,  if  no  reaction  takes 
place,  the  state  remains  constant  and  whenever  the  ith 
reaction  takes  place,  0,-(x)  is  “activated”  and  the  state 
x  is  reset  according  to  (6),  furthermore,  the  probability 
of  the  activation  taking  place  in  an  “infinitesimal”  time 
interval  (t.t  +  dt]  is  A i(x)dt. 


3  Moment  Dynamics 

Given  a  vector  m  =  (m \ ,  m2 , . . . ,  mn)  G  N”  0  of  n  greater 
than  equal  to  zero  integers,  we  define  the  (uncentered) 
moment  of  x  associated  with  m  to  be 


containing  the  top  k  elements  of  p„  which  correspond 
to  the  lower-order  moments  of  interest.  Then,  (11)  can 
be  rewritten  as 


p  =  Ap+Bp  (13) 

where  p  G  RA,  4xoo  denotes  a  matrix  composed 
of  the  first  k  rows  of  the  infinite  identity  matrix 
and  p  G  Rr  contains  all  the  moments  that  appear  in 
the  first  k  elements  of  AocPoo  but  that  do  not  appear  in  p. 

In  this  paper  we  let  the  vector  p  G  RA  contain  all 
the  moments  of  x  of  order  upto  M  G  N>2,  i.e.,  we  con¬ 
sider  an  M.th  order  truncation.  With  this,  the  evolution 
of  vector  p  can  be  written  as  (13)  for  some  matrices 
A  and  B  with  p  G  W  being  a  vector  of  moments  of 
order  M  +  1 .  Our  goal  now  is  to  approximate  ( 1 3)  by  a 
finite-dimensional  nonlinear  ODE  of  the  form 


pim\t)  =  E[x{m\t)],  Vf>  0 


(7)  v=Av  +  B<p(v),  v  =  [vimlKvim2\...ymk)]T  (14) 


where  E  stands  for  the  expected  value  and 


m\  m2 


xv  '  :=  x,  *x. 


(8) 


The  sum  YJ)=\  mj  is  called  the  order  of  the  moment  m. 


Using  results  from  the  SHS  literature,  more  specif¬ 
ically  by  applying  Theorem  1  in  [16]  to  the  SHS 
(3)-(5),  one  can  show  that  the  time  derivative  of  a 
moment  p  is 


A(ra)  =E 


(9) 


where  the  map  (p  :  should  be  chosen  so 

as  to  keep  v(t)  close  to  pit).  This  procedure  is 
commonly  referred  to  as  moment  closure.  We  call 
(14)  the  truncated  moment  dynamics  and  each  element 
<p*rai(jU)  of  <p(p)  the  moment  closure  function  for  the 
corresponding  element  p{m;  in  p. 

The  procedure  we  adopt  to  construct  these  mo¬ 
ment  closure  functions  is  to  first  assume  a  certain 
separable  form  for  each  element  (p^^p)  of  <p(p)  and 
then  matching  time  derivatives  of  p  and  v  at  some 
initial  time  to,  for  every  deterministic  initial  condition 
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of  the  form  x(fo)  =  x  with  probability  one.  The  choice 
of  initial  conditions  is  justified  by  the  fact  that  the 
class  of  deterministic  distributions  forms  a  natural 
basis  for  the  infinite  dimensional  space  containing 
every  possible  state  j,i„  of  ('ll).  Referring  the  reader 
to  [17]  for  further  details  the  above  procedure  leads  to 
the  following  moment  closure  functions.  Assume  that 
each  entry  (pim\/j,)  of  <p  has  the  following  separable 
form 


truncation  M  =  2,  i.e. 

At  =  [JU  ( 1  ’0’°) ,  M  (O’  t  ’0) ,  M  (O’O’  1) ,  M  (2.0,0)  ?  M  (0,2,0)  5  M  (0,0,2)  ? 

At(i-I’0),At(1’0’1),Ai(0’1-1)]r 

(18) 

Table  2  lists  moment  closure  functions  obtained  from 
(16),  which  approximate  different  third  order  moments 
of  x  in  terms  of  the  first  two. 


<P(li,)(v)  =  n(v^)7p  =  vW,  r=(7i (15) 

p= i  v 

with  yp  are  the  unique  solution  of  the  following  system 
of  linear  equations1 


(mp) 
(mi)  ’ 


Vs  = 


(16) 


where  we  define  for  vectors  m  =  (mi , . . . ,  mn)  and  m  = 

{m ,  •  •  • ,  m„) 


r' 1  m  1  . _ r-!n\  r-tr!2  r'ribi 

^(m)  ^mim2  ' 


Then,  when  /t(fo)  =  v(fo),  we  have 


(17) 


d'n(t0) 

dt’ 


d‘v(t0) 
dt ' 


+  ei(x), 


Vi  >  1 


for  every  deterministic  initial  conditions  of  the  form 
x(fo)  =  x  with  probability  one  for  every  integer  x.  In 
the  equation  above,  each  entry  of  the  vector  J  *  is  a 
polynomial  in  x,  whose  degree  exceeds  by  at  least  M 
the  degree  of  the  corresponding  entry  of  the  error  vec¬ 
tor  e,(x).  Thus,  with  increasing  M,  the  truncated  mo¬ 
ment  dynamics  v(f)  provides  a  more  accurate  approxi¬ 
mation  to  the  moments  in  A i(f).  The  striking  feature  of 
the  moment  closure  constructed  is  that  they  are  inde¬ 
pendent  of  the  reaction  parameters  (reaction  rates  and 
stoichiometry)  and  moreover  the  dependence  of  higher- 
order  moment  on  lower  order  ones  is  consistent  with 
x  being  jointly  lognormally  distributed,  in  spite  of  the 
fact  that  the  procedure  used  to  construct  (p  did  not  make 
any  assumption  on  the  distribution  of  the  population. 
For  a  three-specie  reaction  ( n  =  3)  and  a  second  order 


1Ceh  is  defined  as  follows:  Vt,h  6  N>o 


C 


l 

■  — 


l\ 

0, 


l>h 
l  <h 


Table  2:  Moment  closure  function  <p*m'( At)  for  differ¬ 
ent  third  order  moments  Ll i  m  :  with  M  =  2  and  n  =  3. 


At™ 

<p(AHn) 

^ (3,0,0) 

/V2’°’o)Y 

yAi(1’°’°l ) 

^(2,1-0) 

2.oo) \  f  pi1-1’0) \2 

{pio-ho) )  \pbM  ) 

A*(1>1>1) 

pi1- 1,  o)M(o,  1, 1)^(1, 0,1) 

M(l,0,0)w(0,lfl)M(0,0,l) 

4  Modeling  of  Gene  Expression 
with  negative  feedback 

In  this  section  we  consider  negative  feedback  in 
gene  expression  caused  by  the  protein  inhibiting  its 
own  transcription.  Such  auto-regulatory  networks 
are  common  means  of  stabilizing  protein  levels  in 
biochemical  pathways  [18].  This  negative  feedback 
is  realized  by  first  assuming  that  there  is  an  activator 
that  activate  a  gene  to  transcribe  a  mRNA.  Proteins 
are  translated  from  the  gene  at  a  constant  rate.  These 
proteins  can  then  interact  with  the  activator  to  change 
its  shape,  making  it  incapable  of  activation.  Hence, 
more  protein  implies  lesser  number  of  activators,  and 
hence,  lesser  transcription. 

Denoting  the  mRNA,  protein  and  the  activator  by 
Xi,  A2  and  A3,  respectively,  the  above  interactions  can 
be  written  as 

Gene  +  A3  Gene  +  A3  +  Ai , 

Ai  *,  Ai  Aj  +A2,  A2  * 

*-^4A3,  A3^4* 

A3  +  A,  — 1-4  A2  +  * . 


where  i\  denotes  the  factorial  of  £. 
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Hence  the  mRNA  X\  is  transcribed  from  the  gene  at  a 
rate  A,x3,  where  x3  denotes  the  number  of  molecules 
of  the  activator.  The  protein  Xi  is  translated  from  the 
mRNA  at  a  constant  rate  Kp.  Both  mRNA  and  the  pro¬ 
tein  decay  at  rates  dr  and  dp  respectively.  Ka  and  da 
denote  the  birth  and  death  rate  of  the  specie  X3,  respec¬ 
tively.  The  rate  at  which  the  protein  X2  reacts  with  the 
activator  X3  and  takes  it  away  from  the  transcription 
process  is  given  by  k\ .  Note  that  increasing  values  of 
k\  denote  larger  negative  feedback.  If  k\  =  0  then  there 
is  no  negative  feedback.  The  deterministic  chemical 
rate  equations  are  then  given  by 


xiD  =  Krxjn  -  drx iD,  x2d  =  Kpx iD  -  dpX2D  (19a) 
X3 D  =  Ka  -  {da  +  klx2D)x3D  (19b) 


where  xifl,  x2o  and  x3/)  are  continuous  deterministic 
approximates  of  xi,  x2  and  x3,  respectively.  At  steady 
state  we  have 

x3dH.  (20) 


Ka 

dn 


Un 


As  this  deterministic  model  does  not  provide  us  with 
any  information  about  the  stochasticity  in  the  protein 
population,  we  turn  to  a  stochastic  formulation.  The 
above  interactions  can  be  modeled  as  a  SHS  with  reset 
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Xl 
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Xl 
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and  corresponding  transition  intensities  given  by 

Ai(x)  =  Arx3,  X2(x)=drx  1,  A3 (x)  =  Kpx\ ,  (22) 

Mx)  =  dpX2 ,  A5(x)  =  Ka,  Ag  (x)  =  x3da  (23) 
A7(x)  =  kix2x3.  (24) 

Using  (9)  the  time  derivative  of  the  the  vector  /./  which 
is  given  by  ( 1 8)  is 


ft  =  A  +  A/r  +  B/r.  (25) 


for  some  matrices  A,  A,  B  and 


£  =  Lti(0’1’2),  Ai(°’2’1),M(U’1)]3 


Using  Table  2  the  above  system  can  be  closed  using  the 
following  approximations 


(0,1,2) 


/V°,°’2,\  (n{0X1)\2 
y^(0,U0)  J  ^(0,0,1)  J 


A* 


(0,2,1) 


(0,0,2) 

^(0,0,1) 


„  (i,l,0)  „ 

(1,1,1)  =  M 


1,0,0) 


\(H^Y 

J  (n^  J 
(0.1. 1)^(1, 0,1) 
(0,1, 0)^(0, 0,1)  ■ 


(26a) 

(26b) 

(26c) 


The  statistical  variations  in  the  protein  population  can 
now  be  obtained  from  the  closed  system  (25)-(26).  In 
this  paper,  we  use  the  steady  state  coefficient  of  varia¬ 
tion  defined  by 


\ZE[xi(°°)]  -(e[x2(°°)])2 

E[x2(°°)] 

^(0,2,0)  (oo)  _(Al(0, 1,0)  (oo))2 

jU(0,l,  0)(oo) 


(27) 

(28) 


to  quantify  noise  strength  in  the  protein  population.  As 
analytical  solution  of  the  closed  moment  equations  are 
too  complicated  to  be  of  any  use  we  investigate  them 
numerically  by  taking 

Kr  =  1 ,  Kp  =  5  sec~ 1 ,  dr  =  A  sec~ 1 , 
dp  =  .001  sec-1,  da  =  100  sec-1. 


For  comparison  purpose,  we  now  vary  parameters  k\ 
and  Ku  so  as  to  keep  x2d(°°)  =  500  fixed  and  see  its 
effect  on  the  coefficient  of  variation.  This  implies  from 
(20)  that  x3d(°°)  =  .01  and  for  a  given  k\,  Ka  is  always 
chosen  such  that 


Ka  —  da 


1  +  J-X2fl(o°) 


x3dH- 


(29) 


Table  3  lists  the  values  of  CV  for  different  values  of  k\ 
obtained  from  the  steady  states  of  (25)-(26).  Note  that 
CV  is  lower  with  non-zero  values  of  k\  than  with  when 
there  is  no  feedback,  i.e  k\  =  0.  Hence,  we  conclude 
that  such  negative  feedback  reduce  stochastic  varia¬ 
tions  in  the  protein  numbers. 
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Table  3:  Steady  state  coefficient  of  variation  in  the  pro¬ 
tein  population  CV  obtained  from  (25)-(26)  for  differ¬ 
ent  values  of  k\ . 


h 

CV 

0 

.32 

.2 

.28 

.8 

.26 

1.8 

.23 

5  Conclusions  and  Future  work 

An  approximate  stochastic  model  for  chemically 
reacting  systems  was  presented  in  this  paper.  This 
was  done  by  representing  the  population  of  various 
species  involved  in  a  set  of  chemical  reactions  as  the 
continuous  state  of  a  SHS.  With  such  a  representation, 
the  dynamics  of  the  infinite  vector  containing  all 
the  statistical  moments  of  the  continuous  state  are 
governed  by  an  infinite-dimensional  linear  system  of 
ODEs,  which  we  approximate  by  finite-dimensional 
nonlinear  ODEs. 

This  typically  involved  approximating  higher  or¬ 
der  statistical  moments  of  the  population  in  terms  of 
the  lower  order  moments.  Using  derivative  matching 
techniques,  explicit  analytical  formulas  to  construct 
these  approximations  were  provided.  Using  then 
we  showed  that  in  gene  expression,  where  a  protein 
inhibits  its  own  transcription,  the  resulting  negative 
feedback  reduces  stochastic  variations  in  the  protein 
populations. 

An  interesting  line  of  future  work  would  be  look 
at  gene  expression  and  regulation  where  the  mRNA 
transcribed  from  the  gene  is  not  immediately  acces¬ 
sible  for  translation.  This  would  corresponds  to  gene 
expression  in  Eukaryotic  cells  where  the  mRNA  is 
transported  from  the  nucleus  to  the  cytoplasm  where 
translation  takes  place.  This  can  be  incorporated  in  our 
model  by  introducing  another  specie  that  corresponds 
to  the  inactive  mRNA  and  is  being  converted  into  an 
active  mRNA  at  some  constant  rate.  Another  possible 
direction  of  future  work  would  be  to  investigate 
stochasticity  in  gene  cascade  activation  network  where 
protein  expressed  by  one  gene  activates  another  gene 
to  express  a  second  protein  and  negative  feedback 
schemes  within  this  network. 
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